rm(list = ls(all = TRUE))
setwd("C:\\Users\\lenovo\\Desktop\\R")
#6.1 简单线性回归

#加载数据
# 1
## (1)
#示例代码比书上多了两个数据
age <- c(13, 11, 9, 6, 8, 10, 12, 7, 10, 9,
         11, 12, 15, 16, 8, 7, 10, 15,13, 11)
ucr <- c(3.54, 3.01, 3.09, 2.48, 2.56, 3.36, 3.18, 2.65, 3.01,
         2.83, 2.92, 3.09, 3.98, 3.89, 2.21, 2.39, 2.74, 3.36, 3.54, 3.01)
group <- c(rep(0, 8), rep(1, 12))
group <- factor(group, levels = 0:1, labels = c("Normal", "Patient"))
UCR <- data.frame(age, ucr, group)

## (2)
str(UCR)
attr(UCR, "datalabel") <- "UCR in Kaschin-Beck disease children"
attr(UCR, "var.labels") <- c("Age in years", 
                             "Urine creatinine (mmol)", 
                             "Type of children")

## (3)
save(UCR, file = "UCR.rdata")

load("UCR.rdata")
library(epiDisplay)
des(UCR)
summary(UCR)

#数据中无缺失值，可用于直接分析。

#6.1.1 拟合简单线性回归模型
plot(ucr~age, data = UCR, 
       xlab = "Age in years", ylab = "Urine creatine (mmol)")

#建立线性回归模型
mod <- lm(ucr~age, data = UCR)
mod
attributes(mod)
#单独提取某个属性
mod$fitted.values#模型的拟合值


#6.1.2 模型输出结果的解释
summary(mod)
#方差分析
summary(aov(mod))
#总平方和
SST <- sum((UCR$ucr - mean(UCR$ucr))^2)
SST
#残差平方和
SSR <- sum(residuals(mod)^2)
SSR
#回归平方和
SSW <- sum((fitted(mod)-mean(UCR$ucr))^2)
SSW
#决定系数：回归平方和/总平方和
SSW/SST
#调整后的决定系数：加上对变量个数的”惩罚“
#1-(1-SSW/SST)*((n-k)/(n-k-1))
#n为样本量；k为变量的个数
Radj <- 1-(1-SSW/SST)*((18-1)/(18-2))
Radj

#用age的均方（自由度为1）除于残差的均方得到F统计量的值
resid.msq <- sum(residuals(mod)^2)/mod$df.residual 
Fvalue <- SSW/resid.msq
Fvalue
#计算检验变量age效应的p值
#df1为变量age的自由度，df2为残差的自由度
#lower.tail设为FALSE是为了得到曲线下方F值右侧的面积
pf(Fvalue, df1 = 1, df2 = 16, lower.tail = FALSE)

#6.1.3 回归诊断

#在之前的plot上添加回归线
abline(mod)
#加上期望值的点，18号大小，蓝色
points(UCR$age, fitted(mod), pch = 18, col = "blue")
#加上残差
segments(UCR$age, UCR$ucr, UCR$age, fitted(mod), col = "green")
#使用函数residuls提取每个样本的残差值
res <- residuals(mod)
res
sum(res)
sum(res^2)

#好的拟合，残差应当呈正态分布

#定性
#常见的残差正态性检验方法是查看直方图
hist(res)
#对小样本，更好的是作期望的标准正态分值与残差的散点图
#正态qq图
#如散点聚集在一条直线上，就表明残差呈正态分布
qqnorm(res)
qqline(res)

#定量
#Shapiro-Wilk检验
shapiro.test(res)

#作残差与拟合值之间的散点图来看残差的分布模式
#type = "n"表示不显示散点
plot(fitted(mod), res, xlab = "Fitted values", type = "n")
#用text在散点位置处标注了散点的编号
text(fitted(mod), res, labels = rownames(UCR))
#添加蓝色的水平参考线
abline(h = 0, col = "blue")
#看不出明显的分布模式，可认为残差独立于拟合值（期望值）

#另一种方法获取模型的残差诊断图
#直接用pdf~dev.off保存了，右边那个plot有时太小，看不清
par(mfrow = c(2,2))
plot(mod)
par(mfrow = c(1,1))

#6.2 分层线性回归

mod1 <- lm(ucr~age + group, data = UCR)
summary(mod1)
#绘制两条平行的回归曲线
#作散点图，蓝色正常、红色有病
#把因子变量与颜色用ifelse进行转换
col <- ifelse(UCR$group == "Normal", "blue", "red")
col
#设置点的图形，1为空心圆点，19为实心圆点
pch <- ifelse(UCR$group == "Normal", 1, 19)
pch
#画图
plot(ucr~age, data = UCR,
     xlab = "Age in years", ylab ="Urine creatine (nmol)",
     col = col, pch = pch)
#加标签
#“bottomright”, “bottom”, “bottomleft”, “left”, 
#“topleft”, “top”, “topright”, “right”, “center”
legend("topleft", 
       legend = c("Normal children","Diseased children"),
       col = c("blue", "red"), pch = c(1,19))
#提取回归模型的截距、斜率来作回归线
coef(mod1)
#截距相同
b<- coef(mod1)[2]
b
#未患病组
a0 <- coef(mod1)[1]
a0
#患病组
#变量groupPatient取值为1（line13）
#其截距为系数的第一项与第三项之和
a1 <- coef(mod1)[1]+coef(mod1)[3]
a1
#加线
abline(a = a0, b = b, col = "blue")
abline(a = a1, b = b, col = "red")

#在上面，我们假设儿童年龄对尿肌酐含量不受疾病影响
#需要进行统计学检验
#age + group + age:group可以简写为age*group
mod2 <- lm(ucr ~ age + group + age:group, data = UCR)
coef(mod2)

#系数提取
# ucr = coef(mod2)[1]+coef(mod2)[2]*age+
#coef(mod2)[3]*group+coef(mod2)[4]*age*group
a0 <- coef(mod2)[1]
a0
a1 <- coef(mod2)[1] + coef(mod2)[3]
a1
b0 <- coef(mod2)[2]
b0
b1 <- coef(mod2)[2] + coef(mod2)[4]
b1

#绘图
#把因子变量与颜色用ifelse进行转换
col <- ifelse(UCR$group == "Normal", "blue", "red")
col
#设置点的图形，1为空心圆点，19为实心圆点
pch <- ifelse(UCR$group == "Normal", 1, 19)
pch
#画图
plot(ucr~age, data = UCR,
     xlab = "Age in years", ylab ="Urine creatine (nmol)",
     col = col, pch = pch)
#加标签
#“bottomright”, “bottom”, “bottomleft”, “left”, 
#“topleft”, “top”, “topright”, “right”, “center”
legend("topleft", 
       legend = c("Normal children","Diseased children"),
       col = c("blue", "red"), pch = c(1,19))
#加线
abline(a = a0, b = b0, col = "blue")
abline(a = a1, b = b1, col = "red")

#两条回归曲线不平行（斜率不等）
#可能是偶然因素（随机误差）所致
#可能与是否患病有关
#用函数summary查看各个回归系数有无统计学意义
summary(mod2)
#交互作用项age:groupPatient没有统计学意义
#即，这两个斜率不同是偶然因素造成的

#6.3 多重线性回归
#6.3.1 拟合多重线性回归模型
#install.packages("ISwR")
library(ISwR)
data(cystfibr)
?cystfibr
str(cystfibr)
cystfibr$sex <- factor(cystfibr$sex, labels = c("male","female"))
#数据集后5个变量的相关性
cor(cystfibr[,6:10])
#为了简化问题，将fev1设为结果变量
fit1 <- lm(fev1~age+sex+height+weight+bmp, data = cystfibr)
summary(fit1)
#R以哑变量的形式处理分类变量，未被显示的水平（male）为参考组

#6.3.2 多重共线性
#自变量之间由于存在相关性而使模型估计失真或难以准确估计
#比如height、weight与bmp就是典型的有相关性的
cor(cystfibr[,3:5])
#方差膨胀因子VIF
#VIF=1/(1-Rj^2)
#Rj^2为第j个自变量对其余自变量的复相关系数的平方
lm.age <- summary(lm(age ~ sex + height + weight + bmp, 
                     data = cystfibr))
1/(1-lm.age$r.squared)
#VIF>5，存在多重共线性
library(car)
vif(fit1)

#三种解决共线性问题的方法:
#(1)剔除造成共线性的自变量，重新构建模型
#(2)采取主成分回归法，将相关变量整合为一个
#(3)采用逐步回归法，限制有较强相关性的自变量同时进入回归方程

#6.3.3 逐步回归
#drop1可以从包含所有自变量的模型中手动一次剔除一个变量以进行逐步回归
#依据：赤池信息量准则（Akaike Information Criterion， AIC）
#该准则考虑到了模型的拟合优度和用来拟合的参数的数目
#AIC越小，表明模型用较少的参数就获得了足够的拟合优度
drop1(fit1)
#比较AIC，可知应当剔除age来获取最小的AIC
#还可以设置参数，获取似然比检验、F检验的结果
#用F检验比较剔除某个变量后模型与原模型的差异
drop1(fit1, test = "F")

#函数step可以根据AIC值进行逐步回归自动选择“最优”模型，
#参数direction默认为both（向前向后法）
#backword(向后法)，forword（向前法）
fit2 <- step(fit1)
fit2 <- step(fit1, direction = "backward")
fit2 <- step(fit1, direction = "forward")
#可知，只包含变量sex和bmp的模型AIC最小

summary(fit2)

#得到回归系数的置信区间
confint(fit2)

#t检验只表明回归系数与0存在差异
#用F检验来检验变量对模型的统计学意义
anova(fit2)
#结果表明，变量sex、bmp对模型而言都是有意义的
#去掉任何一个变量对模型都会产生显著的影响

#6.3.4 回归诊断
#判断模型是否满足最小二乘法的统计假设
#install.packages("gvlma")
library(gvlma)
gvlma(fit2)

#用epiDisplay汇总模型的主要结果
library(epiDisplay)
regress.display(fit2)

write.csv(regress.display(fit2)$table, file = "mytable.csv")
